% This script is just a simple educational script to demonstrate how to
% access the provided data (Fourier transformed, averaged over the ROIs)
% used for our manuscript
%   --> "Anisotropic Longitudinal Water Proton Relaxation in White Matter 
%        Investigated Ex Vivo in Porcine Spinal Cord with Sample Rotation"

% Author: wallstein@cbs.mpg.de

% NOTE: FOR Room temperature dataset (60° Rot-Angle ... ) --> frequency was
% not correctly adjusted (!)

% Date: 2024-04-18

%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

clear; close all; 

% (1) Load the stored data
load('DataSets_T1_Anisotropy.mat')

% (2) User-Definitions:
str_Temperature         = 'Room';       % 'Room' vs 'Physio'

str_T1Experiments       = 'MP2RAGE';    % 'InvRec' // 'MP2RAGE'
str_IR_RECTorBIR        = 'Rect';       % 'Rect' // 'Bir'

iRotAngle_deg           = 30;       % [0, 10, 20, 30, ... , 80, 90]              

%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% (3) Prepare the (selected) data for plotting

switch str_T1Experiments
    case 'InvRec'
        switch str_IR_RECTorBIR
            case 'Rect'
                DataSets_ROIs       = DataSets_T1_Anisotropy.InversionRecovery.Rect;
            case 'Bir'
                DataSets_ROIs       = DataSets_T1_Anisotropy.InversionRecovery.Bir;
        end
        % + 'Time-Axis' (identical for both types of inversion)
        EvolutionTime_ms    = DataSets_T1_Anisotropy.InversionRecovery.TI_ms;
    case 'MP2RAGE'
        DataSets_ROIs       = DataSets_T1_Anisotropy.MP2Rage;
        EvolutionTime_ms    = DataSets_T1_Anisotropy.MP2Rage.TIs_ms;
end

switch str_Temperature
    
    case 'Physio'
        DataSets_ROIs   = DataSets_ROIs.PhysioTemp;
    case 'Room'
        DataSets_ROIs   = DataSets_ROIs.RoomTemp;
        
end

%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
% (4) Plot the data
% NOTE: Data is stored as '2D-matrix' of dimension 
%           [Rot-Angles x Evolution-Times]

% -> index of rotation angle for plot
idxRotation = find(round(DataSets_T1_Anisotropy.RotationAngles_deg) == iRotAngle_deg);

% NOTE: 90° angle was measured twice --> select which one do you want to
% see here
if(iRotAngle_deg==90)
    % for fist 90° Meas, e.g. idxRotation = idxRotation(1); 
    idxRotation = idxRotation(1);   
end

figure('Name', 'Fig 1: Overview of current dataset')
switch str_T1Experiments
    case 'InvRec'
        plot(EvolutionTime_ms, abs(DataSets_ROIs.IR_curves_ROI_A(idxRotation, :))./max(abs(DataSets_ROIs.IR_curves_ROI_A(idxRotation, :))), ...
            'v', 'Color', 'b', 'DisplayName', 'ROI A'), hold on
        plot(EvolutionTime_ms, abs(DataSets_ROIs.IR_curves_ROI_B(idxRotation, :))./max(abs(DataSets_ROIs.IR_curves_ROI_B(idxRotation, :))),...
            'o', 'Color', 'r', 'DisplayName', 'ROI B'),
        ylim([-0.05 1.05])
        xlabel('Inversion-Time / ms'), ylabel('Normalized Signal')
    case 'MP2RAGE'
        errorbar(EvolutionTime_ms, DataSets_ROIs.MP2RAGE_curves_ROI_A(idxRotation, :), ...
            DataSets_ROIs.MP2RAGE_curves_ROI_A_std(idxRotation, :), ...
            'v', 'Color', 'b', 'DisplayName', 'ROI A', 'MarkerSize', 8,...
            'MarkerFaceColor', 'b'); hold on
        errorbar(EvolutionTime_ms, DataSets_ROIs.MP2RAGE_curves_ROI_B(idxRotation, :),...
            DataSets_ROIs.MP2RAGE_curves_ROI_B_std(idxRotation, :), ...
            'o', 'Color', 'r', 'DisplayName', 'ROI B', 'MarkerSize', 8,...
            'MarkerFaceColor', 'r')
        xlabel('Inversion-Time / ms'), ylabel('Signal / a.u.')
end
title(['Experiment: ' str_T1Experiments ' || Temp: ' str_Temperature])
lgd = legend(); title(lgd, ['Rot-Angle: ' num2str(iRotAngle_deg) '°'])
lgd.Location = 'Best';
set(gca, 'FontWeight', 'bold')


%% (5) Example: IR curve fitting (mono-exponential)


switch str_T1Experiments
    case 'InvRec'
        
        NumOfRotations  = numel(DataSets_T1_Anisotropy.RotationAngles_deg);
        % PreAllocate
        T1_ms_ROIA      = zeros(1,NumOfRotations);
        T1_ms_ROIB      = zeros(1,NumOfRotations);
        
        % Define minimum inversion time for fit:
        TI_minForFit_ms  	= 0;
        [ ~, idx_TImin ]    = find( EvolutionTime_ms >= TI_minForFit_ms, 1, 'first'  );
        
        opt = optimoptions('lsqcurvefit','Algorithm','trust-region-reflective', ...
                'Display','off','MaxIter',1e4,'TolFun',1e-12);
        % Some boundaries    
        lowerBounds     = [0    10      0]; 
        upperBounds     = [2.00 1e4     100];    
        
        % Fit-Function
        Fit_Function    = @(x,xdata)real(x(3)*(1-x(1)*exp(-xdata/x(2))));
        
        for iiRot       = 1:NumOfRotations
            % -----------------------------------------
            % (1): Fit of ROI A
            % Get phase -> last scan (used as reference)
            phaseVal_end            = phase(DataSets_ROIs.IR_curves_ROI_A(iiRot, end));
            Data_forFit_ROI_A       = real(DataSets_ROIs.IR_curves_ROI_A(iiRot,:).*exp(-1i*phaseVal_end));
            
            [coef_FitT1_ROI_A, resnorm_Fit_ROI_A] = lsqcurvefit(Fit_Function, [1.9 700 1.0],...
                EvolutionTime_ms(idx_TImin:end)', (Data_forFit_ROI_A(idx_TImin:end))',...
                lowerBounds, upperBounds, opt);
            % -> save fitresult
            T1_ms_ROIA(iiRot)       = coef_FitT1_ROI_A(2);
            
            % -----------------------------------------
            % (1): Fit of ROI B
            % Get phase -> last scan (used as reference)
            phaseVal_end            = phase(DataSets_ROIs.IR_curves_ROI_B(iiRot, end));
            Data_forFit_ROI_B       = real(DataSets_ROIs.IR_curves_ROI_B(iiRot,:).*exp(-1i*phaseVal_end));
            
            [coef_FitT1_ROI_B, resnorm_Fit_ROI_B] = lsqcurvefit(Fit_Function, [1.9 700 1.0],...
                EvolutionTime_ms(idx_TImin:end)', (Data_forFit_ROI_B(idx_TImin:end))',...
                lowerBounds, upperBounds, opt);
            % -> save fitresult
            T1_ms_ROIB(iiRot)       = coef_FitT1_ROI_B(2);
            
        end
        
        figure('Name', 'Fig 2: Fit IR Experiments')
        plot(DataSets_T1_Anisotropy.RotationAngles_deg, T1_ms_ROIA, 'v', ...
            'Color', 'b', 'DisplayName', 'ROI A'), hold on
        plot(DataSets_T1_Anisotropy.RotationAngles_deg, T1_ms_ROIB, 'o',...
            'Color', 'r', 'DisplayName', 'ROI B')
        xlim([- 10 100]), 
        xlabel('Rotation-Angle / °'), ylabel('T_1 / ms')
        lgd = legend(); title(lgd, ['TI(min) > ' num2str(TI_minForFit_ms) 'ms'])
        lgd.Location = 'Best';
        title(['Experiment: ' str_T1Experiments ' || Temp: ' str_Temperature])
        set(gca, 'FontWeight', 'bold')
end



%% (6) Example: Plot MP2RAGE Datasets


switch str_T1Experiments
    case 'MP2RAGE'
        
        idx_DataSet     = 4;    % NOTE: 13 different datasets were measured 
        kLineShift      = 5;    % arbitrary integer (for demonstration)
        currentData     = DataSets_ROIs.FIDs_AllRotations{idxRotation}{idx_DataSet};
        
        figure()
        % --- first ADC Block
        plot(squeeze(abs(currentData.Fids_TI_1(:,end,currentData.kSpaceCenter))), ...
            'Color', [0.5 0 0], 'DisplayName', ['TI(1) = ' num2str(currentData.TIs_ms(1)) 'ms']), hold on
        plot(squeeze(abs(currentData.Fids_TI_1(:,end,currentData.kSpaceCenter-kLineShift))), ':', ...
            'Color', [0.75 0 0], 'DisplayName', ['k-Line: ' num2str(currentData.kSpaceCenter-kLineShift)]),
        plot(squeeze(abs(currentData.Fids_TI_1(:,end,currentData.kSpaceCenter+5))), '--', ...
            'Color', [1 0 0], 'DisplayName', ['k-Line: ' num2str(currentData.kSpaceCenter+kLineShift)]),
        % --- second ADC Block
        plot(squeeze(abs(currentData.Fids_TI_2(:,end/2,currentData.kSpaceCenter))), ...
            'Color', [0 0 0.5], 'DisplayName', ['TI(2) = ' num2str(currentData.TIs_ms(2)) 'ms']),
        plot(squeeze(abs(currentData.Fids_TI_2(:,end/2,currentData.kSpaceCenter-kLineShift))), ':', ...
            'Color', [0 0 0.75], 'DisplayName', ['k-Line: ' num2str(currentData.kSpaceCenter-kLineShift)]),
        plot(squeeze(abs(currentData.Fids_TI_2(:,end/2,currentData.kSpaceCenter+5))), '--', ...
            'Color', [0 0 1], 'DisplayName', ['k-Line: ' num2str(currentData.kSpaceCenter+kLineShift)]),
        lgd = legend(); lgd.NumColumns = 2; title(lgd, ['Rot-Angle: ' num2str(iRotAngle_deg) '°'])
        lgd.Location = 'Best';
        ylabel('abs(ADC-Signal) / a.u.'), xlabel('#Datapoint of ADC')
        title('Exemplarely MP2RAGE Dataset')
        set(gca, 'FontWeight', 'bold')
end

%% TEST

% figure()
% for iiRot = 1:4
% plot(DataSets_T1_Anisotropy.MP2Rage.PhysioTemp.MP2RAGE_curves_ROI_A(iiRot,:), ...
%     'o', 'DisplayName' ,['Rot: ' num2str(DataSets_T1_Anisotropy.RotationAngles_deg(iiRot)) '°']);
% hold on
% end
% lgd = legend()

        